function varargout = PlotPsychFun_GapsByFreq(varargin)
% PlotPsychFun_GapsByFreq(schedule)
% PlotPsychFun_GapsByFreq(schedule,trialstate)
% PlotPsychFun_GapsByFreq(schedule,trialstate,CurTIDX)
% PlotPsychFun_GapsByFreq(ax,schedule,trialstate,CurTIDX)
% 
% sens = PlotPsychFun_GapsByFreq(...
% [sens,far] = PlotPsychFun_GapsByFreq(...
% 
% Plot psychometric function for gap detection paradigm
% 
% sens is the sens of the psychometric function determined at 50%
% response probability for each frequency.
% 
% far is the false alarm rate
% 
% DJS (c) 2010

if nargin == 1
    schedule = varargin{1};
    ax = gca;
    trialstate = schedule.trialstate;
    CurTIDX = length(trialstate);
elseif nargin == 2
    schedule   = varargin{1};
    trialstate = varargin{2};
    ax = gca;
    CurTIDX = length(trialstate);
elseif nargin == 3
    schedule   = varargin{1};
    trialstate = varargin{2};
    CurTIDX    = varargin{3};
    ax = gca;
elseif nargin == 4
    ax = varargin{1};
    schedule   = varargin{2};
    trialstate = varargin{3};
    CurTIDX    = varargin{4};
end

trialstate = reshape(trialstate,length(trialstate),1);

r = trialstate;

x = strcmpi(schedule.writeparams,'gap_duration');
g = schedule.trials(1:CurTIDX,x);
gaps = unique(schedule.trials(:,x));
if isempty(gaps)
    error('gap_duration parameter not found in schedule.');
end

x = strcmpi(schedule.writeparams,'trial_type');
t = schedule.trials(1:CurTIDX,x);
if isempty(t)
    error('trial_type parameter not found in schedule.');
end

x = strcmpi(schedule.writeparams,'NoiseCF');
f = schedule.trials(1:CurTIDX,x);
freqs = unique(schedule.trials(:,x));
if isempty(freqs)
    error('NoiseCF parameter not found in schedule.');
end

for i = 1:length(gaps)
    for j = 1:length(freqs)
        ms(i,j) = sum(t == 0 & g == gaps(i) & f == freqs(j) & r == -2);  %#ok<AGROW>
        ht(i,j) = sum(t == 0 & g == gaps(i) & f == freqs(j) & r == 1);   %#ok<AGROW>
        fat(i,j) = sum(g == gaps(i) & f == freqs(j) & r == -1);          %#ok<AGROW>
    end
end

ms(find(isnan(ms)))   = 0; %#ok<FNDSB>
ht(find(isnan(ht)))   = 0; %#ok<FNDSB>
fat(find(isnan(fat))) = 0; %#ok<FNDSB>

n1 = ht + ms;
n2 = sum(t == 1);

ht  = ht  ./ n1 * 100;          ht(find(gaps <= 0),:)  = []; %#ok<FNDSB>
fat = fat ./ (n1+fat) * 100;    fat(find(gaps <= 0),:) = []; %#ok<FNDSB>

cr = sum(r == 2) ./ n2 * 100;
fac = sum(ismember(r,[-3,-4])) ./ n2 * 100;

g = gaps(find(gaps > 0)); %#ok<FNDSB>


% fit line and find 50% response for each frequency
sens(:,1) = freqs;
for j = 1:length(freqs)
    i = find(ht(:,j) >= 50,1,'first'); %#ok<FXSET>
    if i == 1, i = 2; end %#ok<FXSET>
    
    [p,s,mu] = polyfit([ht(i-1,j) ht(i,j)],[g(i-1) g(i)],1);
    sens(j,2) = polyval(p,50,[],mu);
end

g = repmat(g,1,size(ht,2));

plot(ax,g,ht,'-o','MarkerSize',10);

hold(ax,'on');

plot(ax,g(:,1),mean(fat,2),':r','LineWidth',2);

x = [min(g(:,1)) max(g(:,1))];

plot(ax,x,[cr cr],':b',x,[fac fac],':r','LineWidth',0.5);

%     plot(ax,[sens sens],[-5 50],'--k',[min(g) sens],[50 50],'--k', ...
%         sens,50,'*k','MarkerSize',10);

hold(ax,'off');

set(ax,'YLim',[-5 105]);
set(ax,'XLim',[min(g(:,1)) max(g(:,1))]);

%     x = sort([g; sens]);
set(ax,'XTick',g(:,1),'XTickLabel',num2str(g(:,1),'%3.1f'));
set(ax,'XScale','log');

xlabel(ax,'Gap Duration (ms)');
ylabel(ax,'% Total');

if nargout >= 1,    varargout{1} = sens;    end
if nargout == 2,    varargout{2} = fac;     end





